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ABSTRACT 

A  model-based  nonlinear  optimal  control  design  is  developed  and  applied  to  a  piezoelectric-driven  nanopo¬ 
sitioning  stage  to  demonstrate  high  accuracy  displacement  control  when  ferroelectric  nonlinearities,  hysteresis 
and  relaxation  behavior  are  present  in  the  piezoelectric  actuator.  The  nonlinear  control  design  is  compared  to 
Proportional-Integral-Derivative  (PID)  control  to  illustrate  performance  enhancements  when  ferroelectric  ma¬ 
terial  behavior  is  directly  incorporated  into  the  control  design.  A  multi-scale  ferroelectric  constitutive  law  is 
used  to  model  the  nonlinear,  hysteretic  and  relaxation  behavior  of  the  piezoelectric  actuator.  The  constitutive 
model  is  implemented  in  a  structural  model  of  a  nanopositioning  stage  that  uses  piezoelectric  stack  actuators. 
Both  moderate  frequency  (100  Hz)  and  quasi-static  (0.2  Hz)  operating  regimes  are  considered  to  address  issues 
associated  with  high  accuracy  atomic  force  microscopy  applications.  The  nonlinear  control  design  includes  an 
open-loop  nonlinear  control  input  determined  from  finite-dimensional  optimal  control  theory  and  linear  pertur¬ 
bation  feedback  determined  from  Proportional  Integral  (PI)  control.  The  hybrid  control  design  reduces  error  by 
over  two  orders  of  magnitude  when  nonlinearities,  hysteresis  and  relaxation  behavior  is  present.  Additionally, 
the  magnitude  of  control  input  is  considerably  smaller  when  the  hybrid  control  design  is  implemented. 

Keywords:  nonlinear  optimal  tracking,  atomic  force  microscopy,  piezoelectric,  relaxation,  perturbation  control 

1.  Introduction 

Piezoelectric  materials  have  been  utilized  in  a  broad  range  of  aerospace,  biomedical,  industrial,  and  automo¬ 
tive  applications.  These  materials  generate  large  forces,  small  displacements  and  fast  response  times  in  response 
to  applied  electric  fields.  These  properties  have  provided  compact  actuators  and  sensors  at  macro-  and  micron- 
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scales.  One  of  the  challenges  in  effectively  utilizing  these  materials  in  devices  relates  to  ferroelectric  nonlin¬ 
earities  and  hysteresis  due  to  reorientation  of  electric  dipoles  at  moderate  to  large  field  levels.  Whereas  this 
behavior  can  be  minimized  by  applying  a  unipolar  field,  material  inhomogeneities  and  residual  fields  often  cause 
minor  loop  hysteresis  from  polarization  switching  at  low  fields.  Ferroelectric  materials  such  as  lead  zirconate 
titanate  (PZT)  are  typically  implemented  in  actuator  applications  due  to  their  relatively  large  piezoelectric 
coupling,  although  moderate  to  large  field  inputs  are  often  necessary  to  compete  with  conventional  actuators 
which  introduces  constitutive  nonlinearities  and  hysteresis.  More  recently,  ferroelectric  relaxor  single  crystals 
(PMN-PT  and  PZN-PT)  have  been  investigated  for  actuator  applications  due  to  their  high  strain  (>  1%)  and 
high  coupling  coefficient  ( k 2  >  90%)  [1-3].  Certain  crystal  cuts  provide  exceptionally  linear  electric  field-strain 
behavior  under  quasi-static  loading  at  low  to  moderate  fields,  although  nonlinearities  and  hysteresis  can  occur 
at  higher  field  levels  and  frequencies.  In  certain  applications,  constitutive  nonlinearities  and  hysteresis  can  be 
compensated  using  classical  PI  or  PID  control,  but  this  can  potentially  lead  to  loss  of  high  accuracy,  band¬ 
width  limitation  and  inefficiencies.  This  is  a  particular  issue  in  nanoscale  position  control  where  piezoelectric 
actuators  are  used  to  accurately  position  a  specimen  for  material  characterization  or  nanofabrication  in  atomic 
force  microscopes  and  scanning  tunneling  microscopes.  Alternatively,  it  is  illustrated  in  [4,  5]  that  the  use  of 
charge-  or  current-controlled  amplifiers  can  essentially  eliminate  hysteresis.  However,  this  mode  of  operation 
can  be  prohibitively  expensive  when  compared  with  more  commonly  employed  voltage-controlled  amplifiers, 
and  current  control  is  ineffective  if  maintaining  DC  offsets  as  is  the  case  when  the  x-stage  of  an  AFM  is  in  a 
fixed  position  while  a  sweep  is  performed  with  the  y-stage.  The  present  analysis  focuses  on  control  development 
for  these  applications  where  nanoscale  position  accuracy  is  paramount  and  charge  control  is  not  feasible. 

The  focus  on  nanoscale  material’s  research  and  device  development  continues  to  increase,  highly  accurate 
nanopositioning  control  becomes  increasingly  important  in  applications  such  image  steering  devices  [6] ,  nanoscale 
material  characterization  [7],  atomic  manipulation  [8],  and  protein  delivery  [9].  Piezoelectric  materials  are  often 
employed  as  actuators  in  nanopositioning  stages  as  well  as  embedded  actuators  within  MEMS  devices  [6].  A 
careful  assessment  of  the  control  design  is  generally  required  to  ensure  adequate  control  is  available  to  meet 
the  performance  objectives  while  accommodating  constitutive  nonlinearities  and  hysteresis.  For  example,  the 
real-time  monitoring  of  biological  processes  using  an  atomic  force  microscope  requires  highly  accurate  position 
control  over  several  microns  which  is  sufficient  to  introduce  constitutive  nonlinearities  and  hysteresis  in  the 
nanopositioning  stage.  In  addition,  low  scan  rates  or  static  displacement  control  can  lead  to  drift  from  relaxation 
mechanisms  thus  reducing  accuracy  in  characterizing  a  specimen  surface  during  quasti-static  measurements. 
Linear  control  designs  are  often  insufficient  to  achieve  the  desired  displacement  accuracy  over  a  broad  range  of 
operating  frequencies  while  compensating  constitutive  nonlinearities  and  hysteresis. 
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An  extensive  amount  of  research  has  focused  on  implementation  of  different  control  strategies  to  compensate 
for  ferroelectric  hysteresis  and  nonlinearities  in  piezoelectric  actuators  -  used  in  nanopositioning  stages  [10-15] . 
Classical  PI  and  PID  designs  have  been  shown  to  be  limited  in  bandwidth  and  are  only  effective  in  small 
scanning  regimes  where  hysteresis  is  minimal.  As  the  nonlinearties  and  hysteresis  increase,  stability  becomes 
an  issue  as  the  control  gains  are  necessarily  increased  to  achieve  the  required  accuracy.  The  application  of 
linear  robust  control  techniques  such  as  TLoo  control  have  provided  improvements  in  bandwidth  and  robustness 
to  nonlinearities  and  hysteresis  [10-12].  However,  when  linear  robust  control  methods  are  implemented,  the 
control  input  focuses  on  compensating  for  the  ferroelectric  constitutive  behavior  which  reduces  the  amount  of 
control  that  can  be  used  to  compensate  for  external  disturbance  loads  and  system  dynamics  of  the  device. 
This  has  led  to  research  focused  on  the  design  and  implemention  of  nonlinear,  model-based  control  designs  to 
accommodate  constitutive  nonlinearities  and  hysteresis  [13,15,16].  In  general,  nonlinear  control  designs  use 
either  a  nonlinear  inverse  filter  to  approximately  linearize  the  material  response  so  that  linear  control  can  be 
applied  or  direct  implementation  of  the  nonlinear  and  hysteretic  constitutive  behavior  within  the  control  design. 
These  two  approaches  are  illustrated  in  Figure  1. 

Previously,  the  nonlinear  control  design  illustrated  in  Figure  1  (ii)  was  implemented  on  a  magnetostrictive  ac¬ 
tuator  to  accurately  track  a  cutting  tool  position  for  machining  out-of-round  automotive  parts  at  high  speed  [16]. 
This  design  utilized  nonlinear  optimal  open-loop  control  to  compensate  for  constitutive  nonlinearities  and  hys¬ 
teresis  while  optimal  perturbation  feedback  was  employed  to  increase  robustness  in  operating  uncertainities.  In 
the  present  analysis,  the  additional  effect  of  relaxation  mechanisms  in  position  tracking  is  considered  and  applied 
to  a  piezoelectric-driven  nanopositioning  stage.  To  illustrate  improvements  provided  by  the  nonlinear  control 
design,  it  is  compared  with  classical  PID  control.  The  nonlinear  control  design  is  shown  to  reduce  tracking  errors 
over  two  orders  of  magnitude  in  operating  regimes  where  nonlinearities,  hysteresis,  and  relaxation  affect  the 
performance.  Additionally,  perturbation  feedback  is  developed  using  PI  control  instead  of  optimal  perturbation 
feedback  previously  employed  [16],  which  provides  a  simplified  control  design  for  real-time  implementation. 

The  development  and  implementation  of  the  nonlinear  control  design  is  presented  as  follows.  In  Section  2  a 
homogenized  energy  model  describing  the  ferroelectric  constitutive  behavior  is  summarized  and  used  to  construct 
partial  differential  equations  (PDE)  and  subsequent  finite  element  model  of  the  piezoelectric  actuator  used  in 
the  nanopositioning  stage.  In  Section  3,  the  control  design  is  developed.  First,  PID  control  is  implemented  to 
illustrate  operating  regimes  where  displacement  accuracy  is  marginal.  The  open  loop  nonlinear  tracking  control 
design  is  subsequently  developed  and  applied  to  the  nanopositioning  stage  to  illustrate  improvements  in  control 
when  nonlinearities,  hysteresis  and  relaxation  behavior  are  incorporated  into  the  control  design.  Lastly,  Section 
4  describes  perturbation  feedback  using  PI  control  to  significantly  improve  robustness  in  operating  uncertainties. 
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(i) 


(ii) 


Figure  1:  (i)  Linear  control  design  employing  an  inverse  filter,  and  (ii)  nonlinear  control  design. 


2.  Piezoelectric  Nanopositioning  Stage  and  Hysteresis  Model 

A  schematic  of  the  nanopositioning  stage  and  atomic  force  microscope  used  in  developing  the  control  design 
is  illustrated  in  Figure  2.  The  side  view  illustrates  the  nanopositioning  stage  and  piezoelectric  actuator  that 
controls  specimen  height  relative  to  the  cantilever.  In  this  mode  of  operation,  the  cantilever  position  is  defined 
by  using  a  photodiode  to  measure  changes  in  a  reflected  laser  beam  and  a  feedback  law  is  used  to  reposition 
the  specimen  to  maintain  constant  surfaces.  A  2-D  scan  in  this  manner  yields  the  surface  topography  of  the 
specimen.  A  top  view  of  the  positioning  stage  is  illustrated  in  Figure  2(ii)  where  two  additional  piezoelectric 
stack  actuators  are  used  to  control  in-plane  displacement.  In  this  design,  in-plane  actuator  coupling  is  typically 
negligible  which  allows  development  of  the  control  design  for  a  single  actuator.  Nanopositioning  stages  may  also 
utilize  piezoelectric  tube  actuators  to  control  in-plane  and  out-of-plane  displacement.  This  design  uses  segmented 
electrodes  where  the  application  of  an  electric  field  forces  the  tube  to  bend  creating  lateral  displacement. 
This  design  has  shown  to  provide  improved  linear  behavior  but  the  non-negligible  coupling  between  in-plane 
displacements  further  complicates  the  control  design  [15].  Whereas  the  control  design  presented  here  is  suitable 
for  accommodating  coupling  present  in  a  piezoelectric-tube  actuator,  the  finite  element  implementation  is  more 
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Figure  2:  Schematic  of  an  atomic  force  microscope  configuration  used  when  constructing  the  nonlinear  control  design, 
(i)  Side  view  of  the  AFM  set-up  with  ^-control  piezoelectric  actuator,  (ii)  Top  view  of  the  positioning  stage  illustrating 
the  configuration  of  the  x  —  y  control  piezoelectric  actuators. 


complex,  e.g.,  see  [17].  The  decoupled  stack  actuator  nanopositioner  design  is  considered  here  to  focus  on  control 
development. 

To  achieve  highly  accurate  displacement  control  in  a  nanopositioning  stage,  it  is  advantageous  to  employ  an 
efficient  and  accurate  ferroelectric  constitutive  law  that  can  be  utilized  in  real-time  control.  Model  development 
employed  in  the  present  analysis  focuses  on  a  homogenized  energy  framework.  The  nonlinear  and  hysteretic 
material  behavior  is  a  manifestation  of  material  relations  at  multiple  length  scales.  At  the  crystal  lattice  level, 
the  applied  electric  field  forces  ions  to  move  and  distort  the  local  lattice  structure  which  can  provide  reversible 
lattice  distortions  or  irreversible  polarization  switching  and  phase  transformations.  At  the  mesoscopic  length 
scale,  regions  of  like  polarization  and  phase  called  domains  form  in  different  orientations  to  reduce  the  total 
energy.  When  a  field  is  applied,  these  domain  structures  reorient  through  a  nucleation  and  growth  process  similar 
to  dislocations  in  metal  plasticity.  This  process  is  strongly  dependent  on  composition  and  defect  structure  and  is 
primarily  responsible  for  the  observed  hysteresis  and  nonlinearities.  At  the  next  larger  length  scale  (grain  level), 
defects  along  grain  boundaries  and  spatial  variations  in  residual  stress  and  local  electric  fields  can  affect  the 
domain  reorientation.  These  effects  yield  the  observed  macroscopic  constitutive  behavior  from  many  grains  fused 
together  to  form  a  ceramic.  In  addition  to  multiple  length  scale  behavior,  ferroelectric  materials  are  dependent 
on  the  drive  frequency  and  temperature.  At  low  frequency  or  static  displacement,  relaxation  mechanisms  can 
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result  in  creep  in  the  nanopositioning  stage.  Rate-dependent  effects  are  incorporated  into  the  constitutive  law 
and  control  design  by  focusing  on  relaxation  behavior  at  low  frequency  and  static  displacement  control  often 
used  in  AFM  scanning  profiles. 

The  multiscale  constitutive  behavior  is  implemented  in  the  control  design  by  utilizing  a  previously  developed 
homogenized  energy-based  model  [18-20].  The  constitutive  model  incorporates  mesocopic  material  behavior 
at  the  domain  or  grain  level  in  a  stochastic  homogenization  framework  to  characterize  macroscopic  material 
behavior.  A  distribution  of  interaction  fields  and  coercive  fields  are  implemented  to  model  the  nucleation  of 
polarization  switching  that  typically  occurs  in  the  presence  of  material  inhomogeneities  and  residual  fields. 
Thermal  relaxation  is  included  in  the  model  by  employing  Boltzmann  relations  to  determine  local  polarization 
switching  when  thermal  energy  affects  the  material  behavior.  Macroscopic  material  behavior  is  determined  by 
homogenizing  the  local  polarization  variants  according  to  the  distribution  of  interaction  and  coercive  fields. 

2.1  Homogenized  Energy  Model 

The  equations  governing  the  homogenized  energy  model  are  summarized  here.  A  detailed  review  of  the  mod¬ 
eling  framework  is  given  in  [19-21].  As  previously  mentioned,  the  control  design  is  developed  for  a  piezoelectric 
stack  actuator  used  in  controlling  the  nanopositioning  stage  for  in-plane  position  control.  The  constitutive  law 
is  therefore  focused  on  uniaxial  loading  of  rod-type  actuators  where  material  coefficients  have  been  reduced  to 
scalar  coefficients  or  distributed  variables  in  the  direction  of  loading.  The  Gibbs  energy  at  the  mesoscopic  length 
scale  is 


G(P,  T)  =  T (P,  T)  —  EP  (1) 

where  lf(P,  T)  is  the  Helmholtz  energy  detailed  in  [18],  T  is  temperature,  E  is  the  electric  field,  and  P  is 
the  polarization.  In  the  one-dimensional  case  considered  here,  the  Helmholtz  energy  function  is  a  double-well 
potential  below  the  Curie  point  Tc  which  gives  rise  to  stable  spontaneous  polarization  with  equal  magnitude  in 
the  positive  and  negative  directions. 

The  effects  of  thermal  relaxation  are  often  present  and  must  be  included  in  the  constitutive  model  to  predict 
creep  in  the  nanopositioning  stage.  This  can  be  accomplished  by  introducing  the  Boltzmann  relation 

fi{G)  =  Ce~GV/kT  (2) 

which  quantifies  the  probability  /i  of  achieving  an  energy  level  G.  Here  kT /V  is  the  relative  thermal  energy 
where  V  is  a  representative  volume  element  at  the  mesoscopic  length  scale,  k  is  Boltzmann’s  constant,  and  the 
constant  C  is  specified  to  ensure  integration  to  unity.  The  inclusion  of  thermal  energy  in  the  energy  function 
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incorporates  the  effect  of  switching  before  minima  in  G  disappear  as  the  temperature  increases.  This  reduces  the 
sharp  transition  of  ferroelectric  switching  as  fields  diametrically  opposed  to  the  polarization  direction  approach 
the  coercive  field. 

The  Boltzmann  relation  gives  rise  to  the  expected  values 

poo  r—Pi 

(P+)  =  /  PpG)dP  ,  <P_>  =  /  Pp(G)dP  (3) 

J  Pi  J  —  oo 

of  the  polarization  associated  with  positively  and  negatively  oriented  dipoles.  Here  ±Pj  are  the  positive  and 
negative  inflection  points  in  the  Helmholtz  energy  definition. 

The  local  polarization  variants  are  defined  by  a  volume  fraction  of  variants  x+  and  x_  having  positive  and 
negative  orientations,  respectively.  The  conservation  relation  a:_  +  x+  =  1  must  hold  for  the  volume  fraction  of 
polarization  variants.  The  kinetic  equations 


X+  =  —p+-X+  +p^  +  x_ 


(4) 


X-  =  -p~  +  X-  +P+-X+ 

govern  the  evolution  of  variants  that  switch  where  the  rate  dependent  behavior  is  determined  by  the  transition 

likelihoods  p _| and  p |_  that  define  the  probabilities  that  polarization  variants  repsectively  switch  into  negative 

or  positive  directions  according  to  the  relaxation  behavior.  The  transition  likelihoods  are  given  by 


1  fp7  e-G^v/kTdP  1  r^+te~G^v/kTdP 

_  x  JPi—e _  _  x  J  —  Pi _  /  r  \ 

P+-  -  T  oo  e_G{EtP)v/kTdp  ,  P-+  -  T(T)  pI+e  G{E,P)V/kTdp  [> 

where  e  is  taken  to  be  a  small  positive  constant.  The  relaxation  time  T  defines  the  time  required  to  switch 
a  polarization  variant.  The  quantity  T2  is  inversely  proportional  to  the  relative  thermal  energy  such  that 
T (T)  =  T\  yJV/kT  where  T\  is  a  time  dependent  material  property.  Additional  details  describing  the  governing 
equations  are  given  in  [19,20].  The  resulting  local  average  polarization  is  quantified  by  the  relation 


P  =  X+(P+)  +  X-(P-).  (6) 

In  cases  where  relative  thermal  effects  are  negligible  (i.e.,  the  relative  thermal  energy  kT /V  is  small)  and 
relaxation  times  are  small  relative  to  the  drive  frequencies,  the  relation  Eq.  (6)  limits  to  the  piecewise  linear 
relation 


P(E  +  Er,  Ec,  0  =  Xe(E  +  E:)  +  SPR 


(7) 
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where  \e  is  the  electric  susceptibility,  Pr  denotes  the  magnitude  of  spontaneous  polarization  and  5  =  1  for 
local  variants  aligned  in  the  positive  direction  and  5  =  —1  for  local  variants  aligned  in  the  negative  direction. 
The  interaction  field  Ej  represents  local  variations  in  the  field  due  to  various  material  inhomogeneities  such 
as  twinned  domain  structures  and  intergranular  residual  stress.  The  spontaneous  polarization  Pr  will  switch 
when  the  diametrically  opposed  effective  field  (Ee  =  E  +  Ej)  reaches  the  coercive  field.  A  semi-infinite  set  of 
polarization  variants,  interaction  fields  and  coercive  fields  Ec  are  used  to  determine  when  each  variant  switches. 

The  macroscopic  polarization  is  computed  from  the  distribution  of  local  variants  from  the  relation 

/OO  /»oo 

/  u(Ec,  Et)  [P  (E  +  En  Ec,  0]  ( t)dEjdEc  (8) 

-oo  J  0 

where  v(Ec ,  Ei)  denotes  the  distribution  of  coercive  and  interaction  fields  and  £  represents  the  initial  distribution 
of  the  local  variants.  The  local  polarization  P  (E  +  Ej\  Ec ,£)  is  determined  from  Eq.  (6)  or  (7)  depending  on 
the  degree  of  relaxation  exhibited  by  the  material.  One  choice  for  the  density  in  which  satisfies  physical  criteria 
is 


v{Ec,Ej)  =  Cie-[ln(Bc/Bc)/2c]2e-i572/2b2  (9) 

where  Ec  is  the  average  coercive  field,  c  quantifies  the  coercive  field  variability,  b  is  the  variance  of  the  interaction 
field,  and  c7  is  a  scaling  parameter.  Techniques  to  identify  general  densities  can  be  found  in  [20].  Numerical 
techniques  to  approximate  Eq.  (8)  and  comparison  to  experimental  results  can  be  found  in  [18]. 

Whereas  the  macroscopic  polarization  is  quantified  by  Eq.  (8),  the  forces  generated  by  the  stack  actuator 
must  be  quantified  for  implementation  within  the  control  design.  This  is  provided  by  the  constitutive  law 

a  =  Ype  +  cDe  -  MP(£)  -  Pr)  ~  h2(P(E)  -  Pr )2  (10) 

representing  uniaxial  stress  in  the  piezoelectric  stack  actuator  where  the  effective  properties  of  the  actuator 
include  Yp  as  the  elastic  modulus  at  constant  polarization,  a  Kelvin-Voigt  damping  parameter  cd ,  £  as  the 
linear  strain  component  in  the  direction  of  loading,  h\  as  the  piezoelectric  coefficient  and  h2  as  the  electrostrictive 
coefficient.  It  is  assumed  that  stress  fields  are  limited  to  the  linear  elastic  regime  where  ferroelastic  switching 
is  negligible.  The  polarization  P(E)  is  computed  using  Eq.  (8)  where  Pr  is  the  initial  macroscopic  remanent 
state  of  the  material.  In  all  simulations  presented,  the  actuator  is  intially  poled  poled  with  a  field  equal  to  2 Ec 
which  results  in  Pr  =  0.206  C/m2. 

Simulations  of  minor  loop  hysteresis  for  the  negligible  relaxation  and  thermal  relaxation  models  are  illustrated 
in  Figure  3  for  the  case  of  zero  applied  stress.  The  negligible  relaxation  model  guarantees  closure  of  minor 


hysteresis  loops  whereas  the  thermal  relaxation  model  illustrates  drift  in  the  polarization  and  strain  behavior. 
A  sinusoidal  waveform  at  1  Hz  was  used  as  the  field  input  in  both  simulations.  The  material  parameters 
associated  with  the  homogenized  energy  model  are  given  in  Table  1.  The  material  parameters  associated  with 
the  density  v(Ec,  Ej)  have  been  adjusted  to  give  approximately  the  same  constitutive  behavior.  This  adjustment 
is  necessary  to  account  for  thermal  energy  in  the  relaxation  model.  The  size  of  the  representative  volume  element 
V  was  selected  to  give  typical  relaxation  behavior  for  an  isothermal  process  and  internal  damping  behavior  at 
room  temperature. 

The  stress  computed  using  Eq.  (10)  includes  linear  stress-strain  behavior  as  well  as  nonlinear  and  hysteretic 
dependence  on  the  electric  field  through  the  P(E)  relation.  It  does  not  include  spatial  dependence.  This  is 
incorporated  in  the  structural  model  in  the  following  section. 

2.2  Structural  Model 

To  facilitate  the  control  design,  the  constitutive  relations  given  by  Eqs.  (8)  and  (10)  are  used  to  develop  a 
system  model  that  quantifies  forces  and  displacements  when  a  electric  held  or  stress  is  applied  to  the  piezoelectric 
stack  actuator  used  in  the  nanopositioning  stage.  The  partial  differential  equation  (PDE)  model  is  first  given 
and  then  formulated  as  set  of  ordinary  differential  equations  (ODEs)  through  a  finite  element  discretization  in 
space.  The  structural  coupling  of  the  nanopositioning  stage  is  modeled  as  a  damped  oscillator  to  account  for 
boundary  conditions  at  the  end  of  the  piezoelectric  actuator.  The  geometry  used  in  constructing  the  structural 
model  is  illustrated  in  Fig.  4. 

A  balance  of  forces  for  the  structural  model  is  given  by  the  relation  [18, 22] 

Table  1:  Parameters  employed  in  the  homogenized  energy  models. 


Thermal  Relaxation 

Rate  Independent 

ci  =  0.54  C/MV2 

ci  =  0.60  C/MV 2 

c=  0.4 

c  =  0.3 

b  =  0.87  MV/m. 

b  =  0.7  MV/m. 

Xe  =  0.008  C/mV 

Xe  =  0.01  C/mV 

T  =  300  K 

V  =  1.25  x  105  nm3 

t  =  100  ns 
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(a) 


(b) 


Figure  3:  Nonlinear  and  hysteretic  homogenized  energy  model  results  for  minor  loop  hysteresis  typically  observed  in 
ferroelectric  materials  under  zero  stress;  comparison  between  the  negligible  relaxation  model  and  thermal  relaxation 
model,  (a)  Electric  field  versus  polarization,  (b)  Electric  field  versus  longitudinal  microstrain. 


d2w  =  dNtot 
dt 2  dx 


(11) 


where  the  density  of  the  actuator  is  denoted  by  p,  the  cross-section  area  is  A  and  the  displacement  is  denoted 
by  w.  The  total  force  Ntot  acting  on  the  actuator  is 


+w 


x=0 


x=L 


Figure  4:  Piezoelectric  stack  actuator  with  damped  oscillator  used  to  quantify  loads  during  scanning  operations  of  an 
AFM  nanopositioning  stage.  Disturbance  forces  along  the  actuator  are  given  by  Fd  and  the  control  input  is  u(t). 
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Table  2:  Model  parameters  for  the  piezoelectric  stack  actuator  and  damped  oscillator. 


Yp  =  60  x  109  N/m2 
hi  =  1.0  x  106  N/C 
h2  =  1.0  x  103  Nm2/C 1 
PR  =  0.3  C/m2 


p  =  7.5  x  103  kg/m3 
A  =  5.07  x  10"4  to2 
cd  =  3.7  x  106  Ns/m 
L  =  0.1  to 


kL  =  3.04  x  106  N/m 
cL  =  3.04  x  102  Ns/m 
m.L  =40  g 


Ntot{t ,  x)  =  YPA ^  +  cDA^-t  +  Fp(E)  +  Fd  (12) 

where  the  elastic  restoring  force  is  given  by  the  first  term  on  the  right  hand  side  of  the  equation  and  Kelvin- Voigt 
damping  is  incorporated  in  the  second  term.  The  linear  elastic  strain  component  in  the  direction  of  loading  is 
defined  by  e  =  '/// .  The  term  Fd  incorporates  external  disturbance  loads  and  the  coupling  force  Fp  represents 
forces  generated  by  an  applied  electric  field  where 

FP(E)  =  A[hi(P(E)  -  Pr)  +  h2{P{E)  -  Prf }  (13) 

and  the  hysteretic  and  nonlinear  E  —  P  relation  is  specified  by  Eq.  (8). 

As  illustrated  in  Fig.  4,  the  boundary  conditions  are  defined  by  a  zero  displacement  at  x  =  0  and  the  balance 
of  forces  at  x  =  L  yields 


dw  d2w 

Ntot(t,L)  =  —kLw(t,  L)  -  cL  —  (t,  L)  -  mL-^-(t,L).  (14) 

The  initial  conditions  are  tu(0,a;)  =  0  and  ///(0,x)  =  0.  Model  parameters  associated  with  the  stack  actuator 
and  damped  oscillator  used  in  the  control  design  are  given  in  Table  2. 

The  strong  form  of  the  PDE  model  given  by  Eq.  (11)  can  be  written  in  the  variational  form  for  finite  element 
implementation.  Multiplication  by  weight  functions  and  integration  by  parts  yields 


,  dw 


Y  A— — I-  cdA 
ox 


d2w 

dxdt 


+  Fp(E )  +  Fd 


kLw(t,  L)  +  cl  (t,  L)  +  mL  ( t ,  L) 


m 


(15) 


where  the  weight  functions  4>{x)  are  required  to  be  differentiable  almost  everywhere  and  satisfy  the  essential 
boundary  condition  w{t1 0)  =  0.  The  weak  form  of  the  model  given  by  Eq.  (15)  is  used  to  obtain  a  matrix  ODE 
system. 
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2.2.1  Approximation  Method 

The  piezoelectric  stack  actuator  model  given  by  Eq.  (15)  is  discretized  in  space  followed  by  a  finite-difference 
approximation.  The  transducer  is  divided  into  N  equal  segments  over  the  length  [0,  L\  with  points  Xi  =  ih ,  i  = 
0,1, ...  ,N  and  step  size  h  =  L/N.  The  spatial  basis  is  comprised  of  the  linear  basis  functions 


1 

h  < 


(x  -  Xi- 1)  , 
{Xi+1  -  x)  , 


Xi- 1  <  X  <  Xi 
Xi  <  X  <  Xi+1 

otherwise 


for  i  =  1, . . . ,  N  —  1.  For  i  =  N,  the  basis  function  is  defined  to  be 


(16) 


0n{x) 


1 

h 


(x  -  XjV-i)  , 
0  , 


XN-l  <  X  <  XN 
otherwise. 


Details  describing  the  basis  functions  can  be  found  in  [18,23]. 

The  approximate  solution  to  Eq.  (15)  is  given  by  a  linear  superposition  of  the  basis  functions 


(17) 


N 

wN{t,  x)  =  VjWMx)  (18) 

0  =  1 

where  Wj  (t)  are  the  nodal  displacement  solutions  along  the  length  of  the  transducer. 

A  second  order  matrix  equation  is  obtained  by  substituting  wN(t,x )  into  Eq.  (15)  with  the  basis  functions 
employed  as  the  weight  functions.  This  gives 


Mw  +  C  £>w  +  Kw  =  FP(E)  b  +  fd  (19) 

where  w (t)  =  . . . ,  lujv(i)],  and  M  G  RNxN ,  Cp  €  RNxN  and  K  G  RNxN  denote  the  mass,  damping  and 

stiffness  matrices.  The  vectors  b  G  Rw  and  fj  G  RN  include  the  integrated  basis  functions  related  to  the  control 
input  and  disturbance  loads,  respectively.  Components  contained  within  the  system  matrices  and  vectors  can 
be  found  in  [16, 18]  and  in  the  Appendix. 

Formulation  of  Eq.  (19)  as  a  first  order  system  yields 

x(f)  =  Ax(f)  +  [B(u)](f)  +  G  (t) 

x(0)  =  x0  (20) 

y(t)  =  Cx(f) 
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where  x(i)  =  [w,w]T  should  not  to  be  confused  with  the  coordinate  x.  The  matrix  A  incorporates  the  mass, 
damping  and  stiffness  matrices  given  in  Eq.  (19)  and  [B(u)](t)  includes  the  nonlinear  input  where  u(t)  is  defined 
as  the  electric  field.  The  initial  conditions  are  defined  by  xo-  The  output  of  the  system  y(t )  is  a  function  of  the 
system  states  according  to  the  matrix  C.  In  the  nanopositioning  stage,  it  is  assumed  that  only  the  displacement 


at  x  =  L  is  observable  which  results  in  C  = 


1  0 


0 


with  dimension  1  x  2 N .  External  disturbances 


are  incorporated  in  G.  The  system  matrix  A,  input  vector  B(w)  and  disturbance  load  G(t)  are 


0  I 

0 

0 

A  = 

-M^K  — M-1C£) 

,  B(u)=Fp(u) 

M-'b 

,  G(i)  = 

1 

£ 

1 

£ 

1 _ 

The  identity  matrix,  with  dimension  N  x  N,  is  denoted  by  I. 

A  temporal  discretization  of  the  system  given  by  Eq.  (20)  is  used  to  numerically  analyze  the  dynamic 
performance.  The  trapezoid  rule  is  adopted  since  it  is  moderately  accurate,  A-stable,  and  requires  minimal 
computer  storage  capacity.  The  trapezoidal  discretization  yields  the  relation  for  each  iteration 


xfc+1  =  Wxfc  +  V  [B(ufc)]  +  [G (4)  +  G(4+r)] 


x(0)  =  x0 


(22) 


where  a  temporal  step  size  At  is  employed  giving  a  discretization  in  time  defined  by  tk  =  kAt.  The  values  X& 
approximate  x(t^).  The  matrices 


W  = 


V  =  At 


-l 


(23) 


are  created  once  for  numerical  implementation  yielding  approximate  solutions  with  0(h2 ,  (At)2)  accuracy.  In 
Eq.  (22),  only  values  for  the  control  input  at  the  present  time  tk  are  included  in  the  temporal  discretization  to 
model  dynamics  of  real-time  feedback  developed  in  the  following  section.  This  approach  is  employed  since  the 
control  input  can  only  be  dependent  on  the  current  and  previous  states. 


3.  Control  Design 

We  first  illustrate  the  performance  and  limitations  of  a  PID  control  design  when  ferroelectric  nonlineari¬ 
ties,  hysteresis  and  relaxation  behavior  are  present  in  the  piezoelectric  actuator.  The  ferroelectric  constitutive 
behavior  is  shown  to  introduce  degradation  in  position  accuracy  when  classical  control  methods  are  employed. 
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Improvements  in  high  accuracy  tracking  is  then  obtained  by  implementing  a  nonlinear  optimal  control  design. 
Actuator  displacements  on  the  order  of  tens  of  microns  are  chosen  for  the  analysis  where  it  has  been  demon¬ 
strated  that  nonlinearities  and  hysteresis  are  significant.  The  reference  tracking  profile  used  in  the  simulations 
is  based  on  a  typical  scan  and  hold  procedure  used  in  nano-  and  micron-scale  materal  characterization.  This 
procedure  requires  that  the  specimen  translate  in  the  axdirection  while  holding  in  the  y-direction.  At  the  end  of 
each  scan,  the  position  is  held  in  the  x-direction  and  increased  or  decreased  a  small  amount  in  the  y-direction 
and  the  procedure  is  repeated  in  the  negative  ^-direction.  Two  scanning  frequencies  are  considered  —  one 
moderate  frequency  (100  Hz)  where  relaxation  is  considered  negligible  and  one  low  frequency  (0.2  Hz)  where  re¬ 
laxation  and  creep  are  significant.  The  constitutive  relations  given  by  Eqs.  (6)  and  (7)  are  employed  to  illustrate 
operating  frequencies  where  hysteresis  and  relaxation  behavior  affect  the  tracking  performance. 

3.1  Proportional-Integral-Derivative  (PID)  Tracking  Control 

We  first  consider  classical  PID  control  design  to  demonstrate  the  need  for  nonlinear  control  when  operating 
in  regimes  where  ferroelectric  hysteresis  and  relaxation  is  significant.  Integral  control  is  first  developed  by 
appending  an  integral  state  to  Eq.  (20)  that  represents  the  error  dynamics  between  the  piezoelectric  actuator 
displacement  and  the  prescribed  displacement  ?’(f). 

The  error  between  the  output  y(t)  and  the  specified  displacement  is  used  to  formulate  the  error  dynamics 
governed  by 


xi(t)  =  Cx(f)  —  r(t). 


(24) 


The  integral  state  augments  the  state  space  system  to  track  the  accumulation  of  error.  This  gives  rise  to  the 
relation 


Xi(t) 

0  C 

Xi(t) 

0 

1 

0 

= 

+ 

u{t)  - 

r(t)  + 

.  *(*) 

0  A 

.  XW 

B 

0 

_  G  (f)  _ 

and  the  control  input  is 


(25) 


u(t) 


K,  K0 


Xi{t) 

x(t) 


(26) 


where  A'/  is  a  scalar  quantity  that  represents  the  integral  control  gain  and  Ko  is  a  vector  with  dimension  1  x  2N 
which  includes  control  gains  on  the  displacements  and  velocities  along  the  length  of  the  piezoelectric  actuator. 
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Therefore  Kq  includes  proportional  control  Kp  on  the  actuator  displacement  w (t)  at  x  =  L  and  derivative 


control  Kp,  on  the  actuator  velocity  w  (t)  at  x  =  L  (i.e.,  Kq  = 


KP  0 


0  KD  0 


0 


)• 


3.1.1  Numerical  Results:  PID  Control  Design 

The  tracking  performance  of  the  nanopositioning  stage  using  PID  control  is  demonstrated  using  scan  rates 
conducted  at  two  different  frequencies.  In  the  first  case,  we  consider  a  moderate  frequency  scan  rate  at  100 
Hz  where  relaxation  behavior  is  expected  to  be  minimal.  The  negligible  relaxation  model  based  on  Eq.  (7) 
is  employed  in  the  constitutive  law  and  the  control  design  for  the  100  Hz  scanning  frequency.  In  the  second 
case,  a  low  frequency  scan  rate  of  0.2  Hz  is  applied  to  the  piezoelectric  actuator  and  the  relaxation  model 
(Eq.  (6))  is  used  to  illustrate  the  additional  effect  of  creep  in  tracking  a  reference  displacement.  In  both 
cases,  the  maximum  displacement  specified  to  the  piezoelectric  actuator  was  40  ym  and  the  disturbance  loads 
G(t)  were  set  to  zero.  Additionally,  it  was  determined  that  a  lumped  paramater  model  (single  finite  element) 
was  sufficient  to  accurately  model  the  piezoelectric  actuator  dynamics  for  a  100  Hz  scan  rate,  thus  the  state 
space  system  reduces  to  scalar  coefficients  with  N  =  1.  The  lumped  parameter  model  assumes  the  a  —  e 
and  E  —  P  are  uniform  along  the  rod  which  has  been  shown  to  be  a  reasonable  approximation  [20] ;  however, 
in  general  this  approximation  is  not  necessarily  valid  for  all  rod  structures  and  operating  regimes.  In  these 
cases,  implementation  of  the  discretized  PDE  is  necessary.  The  control  gains  on  the  integral  state,  actuator 
displacement  and  velocity  at  x  =  L  were  adjusted  to  minimize  the  error  between  the  prescribed  displacement 
and  the  commanded  displacement.  The  values  used  in  the  simulations  were  Kj  =  1.4  x  1012,  Kp  =  lx  108  and 
Kp  =  2000  where  the  units  were  N-m  and  mC-  kV.  Larger  control  gains  resulted  in  no  improvement  in  error 
reduction  and  eventually  led  to  an  unstable  response. 

In  Figure  5(a),  marginal  tracking  performance  is  achieved  using  PID  control  where  the  displacement  am¬ 
plitude  of  40  fxm  resulted  in  a  RMS  error  of  1.2  pm.  Also  note  that  the  maximum  error  occurred  during  the 
scanning  process  and  not  during  static  displacement  control.  This  can  lead  to  erroneous  material  characteri¬ 
zation  as  the  AFM  probe  moves  across  the  specimen  surface.  The  E  —  P  response  associated  with  the  control 
input  is  illustrated  in  Figure  5(b). 

The  tracking  performance  and  constitutive  response  associated  with  the  control  input  at  0.2  Hz  is  illustrated 
in  Figure  6.  When  the  scan  rate  occurs  at  the  lower  frequency,  the  RMS  tracking  error  is  slightly  reduced  to  a 
value  of  0.7  ym.  Whereas  the  error  dynamics  are  similar  to  the  high  frequency  case,  nonzero  error  occurs  under 
static  displacement  control  when  relaxation  behavior  is  present.  PID  control  partially  compensates  for  this 
behavior  as  illustrated  in  Figure  6(b)  where  the  field  input  increases  or  decreases  to  ensure  that  the  polarization 
and  displacement  are  approximately  constant.  Errors  from  relaxation  during  static  displacement  control  can  be 
deleterious  since  the  y— piezo  actuator  (see  Figure  2(b))  must  not  move  during  the  scan  in  the  ^-direction. 
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Figure  5:  Displacement  tracking  performance  using  PID  control  for  a  moderate  frequency  scan  rate  (100  Hz),  (a) 
Tracking  performance  where  the  error  is  defined  as  e(f)  =  C x(t)  —  r(t).  (b)  E  —  P  hysteresis  associated  with  the  control 
input. 


Figure  6:  Displacement  tracking  performance  using  PID  control  at  a  low  frequency  scan  rate  (0.2  Hz),  (a)  Tracking 
performance  where  the  error  is  defined  as  e(f)  =  Cx(f)  —  r(t).  (b)  Polarization  -  field  hysteresis  and  relaxation  behavior 
associated  with  the  control  input. 
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3.2  Nonlinear  Optimal  Tracking 

Development  of  the  nonlinear  optimal  tracking  control  design  follows  a  previous  approach  focused  on  negligi¬ 
ble  relaxation  hysteresis  of  magnetostrictive  actuators  for  vibration  attenuation  of  beam  and  plate  structures  and 
tracking  control  of  rod  structures  [16,24,25].  We  summarize  here  key  equations  associated  with  optimal  tracking 
control  and  its  application  in  compensating  ferroelectric  nonlinearities,  hysteresis  and  relaxation  behavior. 

Optimal  tracking  control  utilizes  a  cost  functional  to  determine  the  optimal  control  input.  The  cost  functional 


J=  \(Cx(tf)-r(tf))TP(Cx(tf)-r(tf))  + [H  -  XT (f)x(f)]  dt  (27) 

penalizes  the  control  input  and  the  error  between  the  piezoelectric  actuator  displacement  and  the  prescibed 
displacement  where  P  penalizes  large  terminal  values  on  the  tracking  error,  H  is  the  Hamiltonian,  and  A (t)  is 
a  set  of  Lagrange  multipliers.  The  Hamiltonian  is 


H  =  ^  [(Cx(t)  -  r(t))TQ(Cx(t)  -  r(t))  +  uT(t)Ru(t)\ 

+\T  [Ax(f)  +  [B(u)](f)  +  G(t)] 


(28) 


where  penalties  on  the  tracking  error  and  the  control  input  are  given  by  the  variables  Q  and  R,  respectively. 

The  minimum  of  the  cost  functional  in  Eq.  (27)  is  determined  under  the  constraint  of  the  differential  equation 
given  by  Eq.  (20).  By  employing  Lagrange  multipliers  an  unconstrained  minimization  problem  is  constructed 
where  the  stationary  condition  for  the  Hamiltonian  yields  the  adjoint  relation  [26, 27] 


A (t)  =  -AtA (t)  -  CTQCx(t)  +  CTQr(t) 


(29) 


and  optimal  control  input 


The  resulting  optimality  system  is 


x{t) 

Ax(t)  +  [B(u)](f)  +  G  (t) 

A  (t) 

-AtA(<)  -  CTQCx(t)  +  CTQr{t) 

x(t0)  =  x0 

A(i/)  =  CTP  (Cx(tf)  —  r(tf)) . 


(30) 


(31) 
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The  force  determined  from  Eq.  (13)  is  included  in  the  input  operator  [B(u)](t)  which  directly  includes  the 
nonlinear  and  hysteretic  E  —  P  relation  as  well  as  relaxation  behavior  within  the  control  formulation.  This 
system  of  equations  results  in  a  two-point  boundary  value  problem  which  presents  challenges  in  obtaining 
a  solution  for  large  systems.  Furthermore,  the  nonlinear  nature  of  the  input  operator  precludes  an  efficient 
Riccati  formulation. 

We  first  not  that  Eq.  (31)  has  the  general  form 


z(f)  =  F(f,  z) 

E0z(f0)  =  [x0,0]T  (32) 

Efz(t/)  =  [0,  0]T 

where  z  =  [x(f),A(f)]T  and 


F(t,z) 


Ax(<)  +  [B(u)](t)  +  G  (t) 

-ATX (t)  -  CrQCx(t)  +  CTQr(t) 


Eq  — 


I 

0 

1 

O 

o 

1 _ 

,  Ef  = 

0 

0 

-CTPC  I 

(33) 


Here  I  denotes  an  identity  matrix  with  dimension  corresponding  to  the  number  of  basis  functions  employed  in 
the  spatial  approximation  of  the  state  variables.  Also  note  that  the  prescribed  displacement  has  been  restricted 


to  r(tf )  =  0. 

The  system  given  by  Eq.  (32)  is  approximated  by  discretizing  the  time  interval  [to  -  tf]  with  a  uniform  step 
size  At  at  the  points  to,  ti ,  ■  ■  •  ,  tjv  =  t/-  The  approximate  values  of  the  state  solutions  and  the  adjoint  at  each 
time  step  are  denoted  by  z0,  •  •  •  ,  zN.  Whereas  several  techniques  are  available  to  approximate  the  solution  to 
Eq.  (32),  such  as  finite  differences  and  nonlinear  multiple  shooting  [28],  a  central  difference  of  the  temporal 
derivative  is  utilized  here  to  yield 


1 

At 


[Zj  +  l  Zj] 


1 

2 


[F  (tj,Zj)  +  F(fy+i,  Zj_|_i)] 


EqZq  —  [xq,  Of 


(34) 


EfZN  =  [0, 0] 

for  j  =  0,  •  •  •  ,  N  —  1.  As  summarized  in  the  Appendix,  the  central  difference  approximation  provides  a  formu¬ 
lation  that  can  be  cast  in  a  analytic  LU  decomposition  to  solve  the  two-point  boundary  value  problem. 
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Equation  (34)  can  be  expressed  as  the  problem  of  finding  Zh  =  [zq,  ■  •  •  ,  z n\  which  solves 


•F(zh)=0.  (35) 

Equation  (35)  includes  the  optimality  system  at  each  time  step  and  the  boundary  conditions.  Details  are  given 
in  [24]. 

A  quasi-Newton  iteration  of  the  form 


where  ^  solves 


7k+i  _  k 
&  u  — 


=  -HO, 


(36) 


(37) 


is  then  used  to  approximate  the  solution  to  the  nonlinear  system  given  by  Eq.  (35).  The  Jacobian  T'{ z fc)  has 
the  form 


H(zh)  = 


where 


So  Rq 

Si  Ri 

Sn-i  Rn-i 

E0  Ef 


I 

0 

1 

A  £B  [<] 

0 

I 

~~  2 

-CtQC  -At 

(38) 


(39) 


for  ,S'j.  The  representation  for  Ri  is  similar. 

Approximate  solutions  using  the  quasi-Newton  iteration  requires  inverting  the  Jabobian  to  obtain  updates  to 
the  states  through  the  iteration  procedure.  When  a  large  number  of  time  steps  and  basis  functions  are  required 
to  represent  the  structure  or  device,  direct  inversion  is  not  feasible.  This  is  avoided  by  implementing  the  analytic 
LU  decomposition.  The  LU  decomposition  is  identical  to  the  formulation  used  in  previous  investigations  [24,25] 
and  is  summarized  in  the  Appendix.  This  gives  rise  to  the  following  representation  of  the  Jacobian 


H(0  =  LU 


(40) 
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where  L  and  U  respectively  denote  lower  and  upper  triangular  matrices.  Updates  to  the  state  and  adjoint 
solution  at  each  time  step  ~^+1  are  obtained  through  direct  solution  of  the  lower  triangular  system  L( 
followed  by  direct  solution  of  the  upper  triangular  system 

3.2.1  Nonlinear  Simulation  Results 

Improvements  in  tracking  control  are  demonstrated  here  by  implementing  the  nonlinear  optimal  control 
design.  The  scanning  frequencies  used  with  PID  control  in  Section  3.1.1  illustrate  the  effectiveness  of  directly 
incorporating  ferroelectric  constitutive  behavior  within  the  nonlinear  control  design.  The  values  used  to  penal¬ 
ized  the  displacement  error  and  control  input  were  Q  =  1  x  1016  and  R  =  1  x  10-8  for  both  the  high  and  low 
frequency  scan  rates.  The  penalty  on  the  final  state  was  P  =  1.4  x  1010. 

In  Figure  7,  the  moderate  frequency  scan  rate  is  simulated  using  the  negligible  relaxation  polarization  kernel 
that  characterizes  typical  nonlinearities  and  hysteresis  at  this  frequency.  The  RMS  error  in  this  case  was 
0.3  /im  for  the  40  /im  displacement  amplitude  illustrating  a  factor  of  four  in  error  reduction  relative  to  PID 
control.  Moreover,  high  performance  tracking  was  achieved  with  a  smaller  applied  field  on  the  piezoelectric 
stack  actuator.  Here  the  maximum  field  input  to  the  stack  actuator  was  less  than  1.5  MV/m  in  comparison  to 
a  maximum  field  input  of  2  MV/m  using  PID  control  at  100  Hz  (Figure  5b). 

In  Figure  8,  the  nonlinear  control  design  is  applied  to  the  piezoelectric  actuator  using  a  scan  rate  of  0.2  Hz. 
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Figure  7:  High  performance  tracking  using  open  loop  nonlinear  optimal  control  with  the  negligible  relaxation  constitutive 
model  at  100  Hz.  (a)  Commanded  displacement  and  error  between  the  actuator  displacement  and  desired  reference  signal, 
(b)  E  —  P  constitutive  behavior  associated  with  the  control  input. 
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Figure  8:  High  performance  tracking  using  nonlinear  optimal  control  to  compensate  for  relaxation  behavior  at  0.2  Hz. 
(a)  Commanded  displacement  and  error  between  the  actuator  response  and  the  reference  signal,  (b)  The  associated 
control  input  and  polarization  relaxation  behavior. 


The  RMS  error  was  0.17  /im  in  this  case.  Whereas  relaxation  behavior  is  present  at  this  operating  frequency 
as  illustrated  in  Figure  8(b),  the  nonlinear  control  design  effectively  compensates  for  this  behavior.  This  is 
clearly  demonstrated  in  regions  of  static  displacement  where  the  control  input  increases  or  decreases  to  ensure 
the  polarization  is  constant  which  gives  rise  to  static  displacement. 

Remark  1:  The  nonlinear  control  problem  relies  on  an  analytic  LU  decomposition  to  invert  the  Jacobian. 
The  Jacobian  is  approximated  by  assuming  the  term  [B(u)]  is  linear  which  provides  constant  values  for 
Si  and  Ri  at  all  time  steps.  This  results  in  suboptimal  control,  but  facilitates  evaluation  of  the  derivative 
of  the  input  operator  with  respect  to  the  adjoint.  Additionally,  numerical  implementation  of  the  analytic  LU 
decomposition  given  in  the  Appendix  requires  careful  selection  of  the  time  step  and  the  penalties  on  the  tracking 
error  and  control  input.  For  example,  increasing  the  number  of  finite  elements  can  create  an  ill-conditioned 
Si  and  Ri  which  can  result  in  numerical  errors  when  inverting  Si.  In  the  quasi-Newton  iteration,  convergent 
solutions  are  obtained  using  a  0.1  ms  time  step  for  the  100  Hz  scan  rate  and  1  ms  time  step  for  the  0.2  Hz 
scan  rate.  A  lumped  parameter  model  was  implemented  to  determine  the  dynamic  performance  for  both  cases. 
This  was  validated  using  PID  control  given  in  Section  3.1  where  several  finite  elements  were  implemented 
and  the  actuator  displacements  were  found  to  be  practically  identical  to  the  lumped  parameter  model.  The 
numerical  approximation  of  the  actuator  velocity,  however,  resulted  in  some  numerical  error  in  the  nonlinear 
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model.  To  achieve  high  performance  tracking,  large  jumps  in  the  actuator  velocity  occur  over  the  scanning 
profile  particularly  at  high  frequency.  This  requires  a  small  step  size  which  is  not  feasible  in  the  quasi-Newton 
method  due  to  numerical  errors  associated  with  ill-conditioned  matrices  S',  and  Ri.  This  issue  is  addressed  in  the 
nonlinear  problem  by  penalizing  actuator  displacement  and  not  the  velocity  in  the  cost  functional.  Additionally, 
the  control  input  computed  during  each  iteration  is  filtered  by  averaging  values  over  two  consectutive  time  steps 
to  improve  convergence  properties.  These  errors  can  be  potentially  reduced  further  by  using  adaptive  step  sizes 
near  regions  of  large  jumps  in  the  velocity,  although  more  work  is  required  to  implement  this  technique. 

Remark  2:  The  negligible  relaxation  homogenized  energy  model  has  been  employed  to  estimate  typical 
ferroelectric  nonlinearities  and  hysteresis  at  the  100  Hz  scan  rate  whereas  the  relaxation  model  was  implemented 
at  low  frequency  to  evaluate  to  effects  of  charge  relaxation  and  creep.  This  approach  has  approximated  the 
actual  rate-dependent  hysteresis  that  may  occur  over  a  broad  frequency  spectrum.  The  relaxation  model  is  valid 
for  isothermal  processes;  however,  thermal  boundary  conditions  and  heat  generation  may  occur  which  requires 
additional  analysis  to  incorporate  these  effects  into  a  multi-scale  framework  over  a  broad  frequency  range.  In 
the  present  analysis,  the  negligible  relaxation  model  and  relaxation  model  are  assumed  to  contain  comparable 
hysteresis  effects  to  illustrate  differences  tracking  performance  when  relaxation  may  be  present. 

3.3  Perturbation  Control 

In  the  previous  section,  a  nonlinear  open  loop  control  design  was  shown  to  provide  accurate  tracking  behavior 
by  compensating  for  nonlinear  and  hysteretic  ferroelectric  material  behavior.  However,  it  is  well  known  that 
such  open  loop  controls  are  not  robust  with  regard  to  operating  uncertainties  such  as  unmodeled  constitutive 
behavior  or  disturbance  loads.  These  issues  are  addressed  by  implementing  perturbation  feedback  determined 
using  PI  control  to  improve  robustness.  The  control  design  presented  here  can  be  extended  to  PID  control; 
however,  additional  control  gains  on  the  actuator  velocity  did  not  improve  tracking  performance. 

The  perturbed  system  consists  of  first-order  variations  in  the  state  space  system  and  initial  conditions 
previously  given  by  Eq.  (20), 


<5x(t)  =  A  <5x(f)  +  B  Su(t)  +  SG(t) 

5y(t)  =  C  Sx(t)  +  w(t)  (41) 

<5x(t0)  =  <5xo, 

where  Su,  <5x  and  Sy  are  first-order  variations  about  the  optimal  input  and  system  states  it*,  x*  and  y* .  The 
nonlinear  input  operator  B(u)  has  been  linearized  around  the  optimal  input  u*(t).  The  perturbation  in  plant 
disturbances  is  denoted  by  SG(t).  In  comparison  to  the  output  in  Eq.  (20),  measurement  noise  w(t)  has  been 


22 


included  in  the  first-order  variation  of  the  output  to  illustrate  robustness  in  the  resulting  design.  Both  the 
variation  in  system  noise  SG(t)  and  the  measurement  noise  w(t)  are  assumed  to  be  random  processes  with  zero 
mean.  The  covariance  for  the  system  and  measurement  noise  is  defined  by  E  [<IG(£)c>Gt(t)]  =  VS(t  —  r)  and 
E  [w(t)wT (t)]  =  W5(t  —  t)  with  the  joint  property  E  [<5G(£)k;t(t)]  =  0. 

Proportional  Integral  (PI)  control  is  implemented  in  the  perturbed  system  in  a  manner  analogous  to  the 
approach  summarized  in  Section  3.1.  An  integrator  state  associated  with  error  in  the  perturbed  states  is 
introduced 


Sxj(t)  =  C  Sx(t)  +  e(t). 

where  e(t)  =  Cx(t)  —  r(f).  This  results  in  an  augmented  perturbation  system  given  by 


(42) 


6xi(t) 

0  c 

6xi(t ) 

0 

1 

0 

= 

+ 

Su(t)  + 

e(t)  + 

Sx(t) 

0  A 

5x(t) 

B 

0 

5G(t) 

The  control  input  is 


Su(t) 


Kj 


K6 


5xj(t) 

5x.(t) 


(43) 


(44) 


where  the  integral  control  gain  Kf  is  a  scalar  coefficient  and  the  proportional  Kp  and  derivative  Kp  control 
gains  are  included  in  the  vector  Kg  similar  to  the  control  gain  presented  in  Section  3.1.  As  previously  mentioned, 
derivative  control  can  be  included,  although  only  PI  perturbation  control  is  considered  here. 

The  nanopositioning  stage  relies  on  position  feedback  using  a  laser  and  a  photodiode.  A  certain  level  of 
sensor  noise  is  expected  to  be  present  which  can  degrade  tracking  performance  from  sensor  feedback.  A  Kalman 
filter  is  used  to  accommodate  measurement  noise  in  the  feedback  sensor,  as  well  as  unmeasured  states,  by 
introducing  the  state  estimator 


8x(t)  =  A Sx(t)  +  B Su(t)  +  L(t)  [ Sy(t )  —  C5x(t)] 
Sx(t0)  =  8x0 


(45) 


into  the  perturbation  system.  Here  (5xo  is  the  estimated  initial  conditions  and  L(£)  is  the  Kalman  gain.  A 
similar  approach  can  be  found  in  [16]. 

Key  equations  used  in  determining  the  Kalman  gain  are  summarized  here;  see  [29,30]  for  additional  details. 
The  Kalman  gain  is  determined  from  the  error  covariance  T(t)  ==  E  [ei(£)e|](r)]  where  the  error  is  associated 
with  first  order  variations  in  the  error  dyanmics 
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ez,(f)  =  Sx(t)  —  5x(t). 


(46) 


By  substituting  Eqs.  (41)  and  (45)  into  Eq.  (46),  the  error  is  determined  by  the  convolution  integral 

eL(t)  =  ®(t,t0)eL(to)  +  f  4>(t,  r )  [SG(t )  —  L(r )w(t )]  dr  (47) 

Jto 

where  <4?(i,  to)  is  the  state  transition  matrix.  A  differential  Riccati  equation  is  obtained  by  manipulating  the 
equation  for  error  covariance;  see  [29,  30]  for  details.  Whereas  this  yields  a  time-dependent  Kalman  gain,  a 
steady  state  solution  is  often  sufficient  to  accurately  estimate  the  system  states.  The  steady-state  Kalman  gain 
is  employed  here  which  is  determined  from  the  algebraic  Riccati  equation 

AT(t)  +  T (t)AT  -  T(t)CTW~1CT(t)  +  V  =  0  (48) 

with  the  initial  conditions  Y(fo)  =  To-  The  steady-state  Kalman  gain  is  then  given  by 

L=TCtW~1.  (49) 

In  summary,  the  matrix  system 


Sxi  (■ t ) 

0 

0  C 

5xi{i) 

e(t) 

Sx(t) 

= 

-BRA 

A  -BK0 

5x(t) 

+ 

6G{t) 

5x(t) 

-BRA 

LC  A  -  BK0  LC 

Sx(t) 

Ldw(<) 

(50) 

Sxi(t0)  =  0 
Sx(t0)  =  Sx0 

Sx(t0)  =  Sx0. 

governing  the  dynamics  of  the  perturbed  state  space  system  with  PI  control  and  Kalman  filter  is  solved  to 
obtain  the  system  response  when  perturbation  feedback  of  the  estimated  states  is  implemented. 

3.3.1  Numerical  Results 

The  hybrid  nonlinear  control  method  using  PI  perturbation  feedback  is  implemented  numerically  using  a 
temporal  discretization  method  similar  to  Eq.  (22).  Since  the  error  e(t)  is  only  known  at  the  present  and  past 
time  steps,  the  temporal  discretization  of  Eq.  (50)  is  formulated  as 
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(51) 


SXk+1  =  W 5Xk  +  VEk  +  i  V  [Nfc  +  Nfc+1] 

A-(O)  =  A0 

r  i  t  r  -|  t  r  j 

where  <5Afe  =  Sxik  6xk  6xfc  ,  E k  =  ek  0  0  and  Nk  =  0  SGk  wk  for  each  time  step 

tk.  The  matrices  W  and  V  are  similar  to  those  used  in  Eq.  (22)  except  the  A  matrix  contains  the  input 

vector,  control  gains,  and  Kalman  gains  given  in  the  form  of  the  matrix  in  Eq.  (50).  This  form  of  the  temporal 

discretization  gives  a  direct  comparison  of  performance  improvements  relative  to  the  PID  control  given  in 

Section  3.1.  The  control  gains  were  Kj  =  5  x  1017  and  Kp  —  1  x  1011  using  the  units  N-m  and  mC-kV.  The 

values  used  in  determining  the  Kalman  gains  were  V  =  5  x  106  and  W  =  5  x  10~6  while  the  measurement  noise 

was  a  random  signal  with  amplitude  1  /uni.  The  first  order  variation  in  disturbance  noise  was  set  to  zero. 

Performance  enhancement  using  perturbation  feedback  is  illustrated  by  considering  unmodeled  constutitive 
nonlinearities  and  hysteresis.  This  is  conducted  by  modifying  constitutive  parameters  previously  given  in  Table  1. 
Only  the  negligible  relaxation  model  is  considered;  however,  similar  tracking  behavior  was  obtained  with  the 
thermal  relaxation  model.  The  modifications  in  the  parameters  resulted  in  deviations  in  the  constitutive  behavior 
as  illustrated  in  Figure  9(a)  where  the  control  input  previously  computed  from  the  nonlinear  open  loop  control 
design  was  simulated.  These  variations  in  the  constitutive  behavior  have  a  significant  effect  on  tracking  control 
thus  motivating  the  need  for  perturbation  feedback;  e.g.,  see  Figure  9(b).  The  RMS  error  in  this  case  was  2  /um. 

Perturbation  feedback  provides  considerable  improvements  in  tracking  performance  relative  to  nonlinear 
open  loop  control  when  unmodel  constitutive  behavior  is  present;  e.g.,  see  Figure  10.  The  RMS  error  is  reduced 
to  2  nm  for  the  40  /um  displacement  which  provides  a  factor  of  600  improvement  in  error  reduction  relative  to 
PID  control.  The  open  loop  control  input  u*  and  the  perturbation  control  input  Su*  are  also  given  to  illustrate 
the  relative  magnitudes  of  the  field  input.  Based  on  the  simulations,  the  perturbation  input  is  relatively  small 
which  verifies  that  linearizing  the  input  operator  in  Eq.  (41)  is  a  reasonable  approach. 

Remark  3:  The  perturbation  feedback  control  has  focused  on  classical  PI  control  which  provides  a  simplified 
control  design  for  experimental  implementation.  Derivative  control  gains  did  not  improve  tracking  performance 
which  is  believed  to  be  due  to  small  numerical  errors  in  computing  actuator  velocity.  As  previously  mentioned, 
large  jumps  in  velocity  particularly  at  high  frequency  require  modifications  to  the  nonlinear  model  to  im¬ 
prove  numerical  accuracy  in  computing  actuator  velocity.  Although  non-zero  derivative  control  gain  introduced 
additional  tracking  errors  in  the  present  model,  the  control  design  using  PI  perturbation  feedback  provides 
exceptional  improvement  in  tracking  control  over  classical  PID  control. 
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Figure  9:  Evaluation  of  degradation  in  tracking  control  when  unmodeled  constitutive  behavior  is  present,  (a)  The  original 
parameters  used  in  the  constitutive  model  are  modified  and  the  E—  P  behavior  is  shown,  (b)  Open  loop  nonlinear  control 
illustrates  large  tracking  errors. 


Figure  10:  The  hybrid  control  design  using  PI  perturbation  control  illustrates  significant  improvements  in  reducing 
tracking  errors,  (a)  The  optimal  open  loop  control  u*  and  the  additional  perturbation  feedback  8u* .  (b)  High  performance 
tracking  behavior  using  perturbation  feedback.  Note  that  the  error  scale  has  been  significantly  reduced  relative  to 
Figure  9. 
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4.  Concluding  Remarks 

A  hybrid  nonlinear  optimal  control  design  has  been  applied  to  a  nanopositioning  stage  for  both  quasi-static 
and  high  frequency  (100  Hz)  scanning  rates  for  applications  where  nanoscale  resolution  is  critical.  Ferroelectric 
material  behavior  was  incorporated  into  the  control  design  using  a  homogenized  energy  model  that  accounts  for 
local  residual  interaction  fields  and  variations  in  the  coercive  field.  Material  parameters  were  chosen  to  represent 
typical  hysteresis  behavior  for  operating  frequencies  that  range  from  tens  to  hundreds  of  Hertz.  Relaxation 
behavior  was  incorporated  into  the  constitutive  model  to  account  for  creep  during  quasi-static  actuation.  The 
hybrid  control  design  significantly  improved  reduction  in  tracking  errors  relative  to  PID  control  where  the  error 
was  reduced  by  over  two  orders  of  magnitude  for  the  high  frequency  scan  rate.  Furthermore,  the  improvement 
in  tracking  performance  was  achieved  using  a  smaller  field  input  which  may  have  implications  in  applications 
were  energy  efficiency  is  critical. 
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Appendix 

A.  Mathematical  Relations 

Matrix  and  Vector  Relations 

The  matrices  and  vectors  given  in  Section  2.2  are  formulated  using  linear  basis  functions.  The  approximate 
displacements  in  Eq.  (18)  are  substituted  into  the  variational  form  of  the  dynamic  equation  Eq.  (15)  which  gives 

N  ,.L  N  JV  ~L  fL 

y^u)j(f)  /  pA(j)j(j)jdx  +  wj ( t )  /  cpAfi'^dx  +  ^  Wj(t)  /  YM Acjj'^dx  =  Fmag(H)  /  dx 

j=l  “'O  j  =  i  A)  j=1  Jo  Jo 

+  /  Fd(j>idx  -  (kLwN(t)<j>N(L)  +  cLwN(t)<j>N(L)  +  mLwN(j>N(L ))  4>N(L) . 

Jo 

The  global  mass,  stiffness,  and  damping  matrices  have  the  components 


/  pA<j>i(j)jdx  , 

Jo 

pA(j>i<f>jdx  +  m,L  , 


i  ^  N  or  j  N 
i  =  N  or  j  =  N 
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[CD\ij 


coA^tf/Ax 


J  CjjAcp'^'jdx  +  cL  , 


i  ^  N  or  j  ^  N 


i  =  N  or  j  =  N 


L 

Ym  Atfrf'jdx  ,  i^N  or  j^N 

L 

YM Atp'ifijdx  +  k,L  ,  i  =  N  or  j  =  N 


whereas  the  force  vectors  have  the  components 


[  fiidx. 
o 

Details  regarding  the  formulation  of  these  matrices  can  be  found  in  Chapter  8  of  [18]. 

LU  Decomposition 

The  Jacobian  given  by  Eq.  (38)  is  written  in  terms  of  a  LU  decomposition  to  avoid  matrix  inversion  when 
considering  large  systems.  The  analytic  form  of  the  decomposition  used  in  the  simulations  is 

So 

Si 


[fd]j  =  /  Fd<t>idx  ,  [b];  = 

Jo 


L  = 


Sn-  i  0 

N- 2  N—l 
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This  form  of  the  matrix  system  reduces  the  computation  to  solving  2N  x  2 N  systems  for  each  time  step 


separately. 
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